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ABSTRACT 

We compare the properties of ordinary strong magnetohydrodynamic (MHD) turbulence in a strongly 
magnetized medium with the recently discovered viscosity-damped regime. We focus on energy spectra, 
anisotropy, and intermittency. Our most surprising conclusion is that in ordinary strong MHD turbu- 
lence the velocity and magnetic fields show different high-order structure function scalings. Moreover 
this scaling depends on whether the intermittency is viewed in a global or local system of reference. 
This reconciles seemingly contradictory earlier results. On the other hand, the intermittency scaling 
for viscosity-damped turbulence is very different, and difficult to understand in terms of the usual phe- 
nomenological models for intermittency in turbulence. Our remaining results are in reasonable agreement 
with expectations. First, we find that our high resolution simulations for ordinary MHD turbulence show 
that the energy spectra are compatible with a Kolmogorov spectrum, while viscosity-damped turbulence 
shows a shallow k~ x spectrum for the magnetic fluctuations. Second, a new numerical technique con- 
firms that ordinary MHD turbulence exhibits Goldreich-Sridhar type anisotropy, while viscosity-damped 
MHD turbulence shows extremely anisotropic eddy structures. Finally, we show that many properties 
of incompressible turbulence for both the ordinary and viscosity-damped regimes carry over to the case 
of compressible turbulence. 

Subject headings: turbulence — MHD — ISM: general - ISM: structure 



1. INTRODUCTION 

The interstellar medium (ISM) shows density and veloc- 
ity statistics that indicate the existence of strong turbu- 
lence (Armstrong, Rickett, & Spangler 1995; Lazarian & 
Pogosyan 2000; Stanimirovic & Lazarian 2001; see also the 
review by Cho, Lazarian, & Vishniac 2003 and references 
therein). Understanding interstellar turbulence is essen- 
tial for many astrophysical processes including star forma- 
tion, and cosmic ray transport (see reviews by Vazquez- 
Semadeni et al. 2000; Lazarian, Cho, Yan 2002 and refer- 
ences therein). Moreover, the importance of magnetohy- 
drodynamic (MHD) turbulence is not limited to interstel- 
lar processes. For instance, properties of MHD turbulence 
may be essential for describing gamma ray bursts (see 
Lazarian et al. 2003) and magnetic reconnection (Lazarian 
& Vishniac 1999). 

Kolmogorov theory (1941) is the simplest model for in- 
compressible hydrodynamic turbulence. The main pre- 
diction of the theory is the velocity-separation relation 
vi ~ Z 1 / 3 , or energy spectrum E(k) ~ fc~ 5 / 3 , where I is 
the separation between two points (or, we regard it as 
eddy size) and k ~ l/l is the wave number. In general, 
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Astrophysical turbulence differs from its lab- 



oratory counterpart in that it is usually magnetized, com- 
pressible, and the gas can be partially ionized. However, 
the single most important difference is that typical large 
eddy scales in MHD turbulence are usually many orders 
of magnitude larger than any relevant microphysical scale. 
This allows us to sidestep many of the issues that make 



laboratory plasma physics so complicated. 

Including magnetic field effects has a dramatic effect on 
strong turbulence. Our current understanding of MHD 
turbulence owes a great deal to fundamental results ob- 
tained by Shebalin, Matthaeus, & Montgomery (1983), 
Higdon (1984), Matthaeus et al. (1998) and other re- 
searchers (see a more complete list in Cho, Lazarian & 
Vishniac 2003). Recently Goldreich & Sridhar (1995; here- 
after GS95) proposed a model for incompressible MHD 
turbulence, which was later supported by numerical simu- 
lations (Cho & Vishniac 2000; Maron & Goldreich 2001). 
The GS95 model predicts a Kolmogorov energy spectrum, 
E(k) ~ fc~ 5 / 3 , for both velocity and magnetic fields. The 
major difference between MHD and hydrodynamic turbu- 
lence is the anisotropy of the eddy structures - eddies are 
statistically isotropic in hydrodynamic turbulence while 
eddies show scale-dependent anisotropy (i.e. smaller ed- 
dies are more anisotropic) in the MHD case (see Cho, 
Lazarian, & Vishniac 2003). 

How useful is studying incompressible turbulence? Does 
it have anything to do with turbulence in realistic com- 
pressible flows? Compressible MHD turbulence has been 
studied intensively by different researchers (for a review 
see Cho & Lazarian 2003c and references therein). A dis- 
cussion of the extent to which the incompressible GS95 
scaling is applicable to realistic compressible fluids may 
be found in the original GS95 paper. More recently Lith- 
wick & Goldreich (2001) extended GS95 model for high-/3 
(/? = Pgas/Pmagnetic) plasmas, i.e. for the case when gas 
pressure is larger than magnetic pressure. They also made 
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a conjecture about the scaling of slow modes in low-/? plas- 
mas. Cho & Lazarian (2002) studied turbulence in low-/? 
plasmas and obtained numerical scaling relations for both 
compressible and Alfvenic modes. Finally, both the high- 
and low-/3 plasma cases are covered in Cho & Lazarian 
(2003a, b). They showed that the Alfvenic and slow modes 
follow the GS95 scaling not only for Mach numbers much 
less than unity, but also for much larger Mach numbers, 
as long as the motions were still sub- Alfvenic. Fast modes 
were shown to be decoupled from the rest of the modes 1 
and isotropic. The work on compressible turbulence shows 
that incompressible simulations arc meaningful in terms of 
representing most features of realistic flows. Consequently 
in this paper we concentrate on incompressible turbulence. 

In this paper, we focus on the effects of viscosity. In 
strong hydrodynamic turbulence energy is injected at a 
scale L, and cascades down to smaller scales without sig- 
nificant viscous losses until it reaches the viscous damping 
scale Idv The Kolmogorov energy spectrum applies to the 
inertial range, i.e. all scales between L and ldv This sim- 
ple picture becomes more complicated when we deal with 
MHD turbulence because there are two dissipation scales - 
the velocity damping scale ld v and the magnetic diffusion 
scale Idm, where magnetic structures are dissipated. In 
fully ionized collisionless plasmas (e.g. the hottest phases 
of the ISM), ldv is less than an order of magnitude larger 
than Idm, but both scales are very small. However, in par- 
tially ionized plasmas (e.g. the warm or cold neutral phase 
of the ISM), the two dissipation scales are very different 
and ldv S> Idm- In the Cold Neutral Medium (see Draine 
& Lazarian 1999 for a list of the idealized phases) neu- 
tral particle transport leads to viscous damping on a scale 
which is a fraction of a parsec. In contrast, in these same 
phases ld m ~ lOOfcm. 

This has a dramatic effect on the energy cascade model 
in a partially ionized medium. When the energy reaches 
the viscous damping scale ldv , kinetic energy will dissipate 
there, but the magnetic energy will not. In the presence of 
dynamically important magnetic field, Cho, Lazarian, & 
Vishniac (2002b; hereafter CLV02b) reported a completely 
new regime of turbulence below the scale at which viscosity 
damps kinetic motions of fluids 2 . They showed that mag- 
netic fluctuations extend below the viscous damping scale 
and form a shallow spectrum Eb(k) <~ fc _1 . The spectrum 
is similar to that of the viscous-convective range of a pas- 
sive scalar in hydrodynamic turbulence (see, for example, 
Lesieur 1990). In addition they showed that turbulence 
in the new regime is very anisotropic and intermittent. 
Here, we compare ordinary and viscosity-damped MHD 
turbulence. We mainly focus on the energy spectra, eddy 
anisotropy, scale-dependent intcrmittency, and intcrmit- 
tency from high order velocity statistics. 

In what follows, we briefly consider a theoretical model 
for viscosity-damped turbulence (§2). A detailed theo- 
retical study is given in Lazarian, Vishniac & Cho (2003; 
hereafter LVC03). We discuss our numerical method in §3. 
We present our results on anisotropies and intermittency 
in §4, while results on higher-order statistics are given in 



§5. A discussion and the summary are presented in §6 and 
§7 respectively. 

2. A THEORETICAL MODEL FOR VISCOSITY-DAMPED 
MHD TURBULENCE 

Following the usual treatment of ordinary strong MHD 
turbulence, we define the wavenumbers fcii and k± as the 
components of the wavevector measured along the local 
mean magnetic field and perpendicular to it, respectively. 
In this Here the local mean magnetic field is the direction 
of the locally averaged magnetic field, which depends not 
only on the location but also the volume over which the 
average is taken. See Cho & Vishniac (2000) and Cho, 
Lazarian, & Vishniac (2002a; hereafter CLV02a) for de- 
tails. 

Lazarian, Vishniac, & Cho (LVC03) proposed a theo- 
retical model for viscosity-damped MHD turbulence. We 
summarize the model as follows. 

Since there is no significant velocity fluctuation below 
Idv, the time scale for the energy cascade below ld v is fixed 
at the viscous damping scale. Consequently the energy 
cascade time scale t cas is scale-independent below ld v and 
the requirement for a scale independent energy rate rate 
yields 

(1) 



bf /teas 



bi <~ constant, or Ef,(k) <~ k , 

where kE b (k) ~ bf. 

In LVC03, we assume that the curvature of the magnetic 
field lines changes slowly, if at all, in the cascade: 



cm <~ constant. (2) 

This is consistent with a picture in which the cascade is 
driven by repeated shearing at the same large scale. It 
is also consistent with the numerical work described in 
CLV02b, which yielded a constant ku throughout the vis- 
cously damped nonlinear cascade. A corollary is that the 
wavevector component in the direction of the perturbed 
field is also approximately constant, so that the increase 
in k is entirely in the third direction. 

The kinetic spectrum depends on the scaling of inter- 
mittency. In LVC03, we define a filling factor </>;, which is 
the fraction of the volume containing strong magnetic field 
perturbations with a scale I ~ k' 1 . We denote the veloc- 
ity and perturbed magnetic field inside these subvolumes 
with a "~" so that 



<f>lVl, 



(3) 



and 



bf 



W- (4) 
We can balance viscous and magnetic tension forces to find 

max[6;fc c , B k^ c }bi ~ k c b x , (5) 
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vi 



where k r 



1 /ldv and fc|i _ c is the parallel component 
of the wave vector corresponding to the perpendicular 
component k c . We used the Goldreich-Sridhar scaling 
{Bok\\ c ~ bi k c ) and bi > bi to evaluate the two terms 
in the square braces. Motions on scales smaller than 

Prior to this research it was erroneously believed that the linear modes are strongly coupled in the case of MHD turbulence and this was 
claimed to be the reason for a rapid decay of turbulent motions. A recent study by Vestuto, Ostriker & Stone (2003) confirmed our findings 
on a marginal coupling of compressible and incompressible MHD motions. 

2 Further research showed that there is a smooth connection between this regime and small scale turbulent dynamo in high Prandtl number 
fluids (see Schckochihin et al. 2002). 
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will be continuously sheared at a rate t~ . These struc- 
tures will reach a dynamic equilibrium if they generate a 
comparable shear, that is 

y <~ t^ 1 ~ constant. (6) 

Combining this with equation (5), we get 

<t>l ~ Kl (7) 

and 

E v (k) ~ k-\ (8) 

Note that equation (5) implies that kinetic spectrum would 
be E v (k) <~ k~ 5 if </>;=constant. 

3. METHOD 

3.1. Numerical Method 

We have calculated the time evolution of incompressible 
magnetic turbulence subject to a random driving force per 
unit mass. We have adopted a pseudo-spectral code to 
solve the incompressible MHD equations in a periodic box 
of size 2n: 

^ = (V x v) x v - (V x B) x B + z/V 2 v + f + VP', (9) 

— =Vx (vxB)+r/V 2 B, (10) 

V-v = V-B = 0, (11) 

where f is a random driving force, P' = P/p + v- v/2, v is 
the velocity, and B is magnetic field divided by {Airp) 1 / 2 . 
In this representation, v can be viewed as the velocity mea- 
sured in units of the r.m.s. velocity of the system and B as 
the Alfven speed in the same units. The time t is in units 
of the large eddy turnover time (<~ L/V) and the length 
in units of L, the scale of the energy injection. In this sys- 
tem of units, the viscosity v and magnetic diffusivity r\ are 
the inverse of the kinetic and magnetic Reynolds numbers 
respectively. The magnetic field consists of the uniform 
background field and a fluctuating field: B = B + b. 
We use 21 forcing components with 2 < k < Vl2, where 
wavenumber k is in units of L _1 . Each forcing component 
has correlation time of one. The peak of energy injection 
occurs at k w 2.5. The amplitudes of the forcing compo- 
nents are tuned to ensure v fa 1 We use exactly the same 
forcing terms for all simulations. The Alfven velocity of 
the background field, B , is set to 1. In pseudo spectral 
methods, the temporal evolution of equations (9) and (10) 
are followed in Fourier space. To obtain the Fourier com- 
ponents of nonlinear terms, we first calculate them in real 
space, and transform back into Fourier space. The av- 
erage helicity in these simulations is not zero. However, 
previous tests have shown that our results are insensitive 
to the value of the kinetic helicity. In incompressible fluid 
dynamics P' is not an independent variable. We use an 
appropriate projection operator to calculate VP' term in 
Fourier space and to enforce the divergence-free condition 
(V • v = V • B = 0). We use up to 256 3 collocation 
points. We use an integration factor technique for kinetic 
and magnetic dissipation terms and a leap-frog method 
for nonlinear terms. We eliminate the 2At oscillation of 
the leap-frog method by using an appropriate average. At 
t = 0, the magnetic field has only its uniform component 



and the velocity field is restricted to the range 2 < k < 4 
in wavevector space. 

For the ordinary turbulence, we mostly use 8-th order 
hyper-viscosity and hyper-diffusion, so that the viscosity 
and diffusion terms in the equations (9) and (10) become 

-^ 8 (V 2 ) 8 v and -r? 8 (V 2 ) 8 B, (12) 

respectively. Here, i/g (= rj§) is adjusted in such a way that 
the dissipation cutoff occurs right before k ~ N/3, where 
N is the number of grids in each spatial direction. This 
way, we can avoid the aliasing error of pseudo-spectral 
method. 

For viscosity-damped turbulence, we mostly use a phys- 
ical viscosity (y = 0.015) and a third order hyper-diffusion 
for magnetic field, so that the dissipations in the equation 
(10) is replaced with 

-%(V 2 ) 3 B. (13) 

However, for the run 256PP-B 1, we use a physical diffu- 
sion for magnetic field. 

We list parameters used for the simulations in Table 1. 
The run 384PH3-Bol is exactly the same as the run with 
the same name in CLV02b. We use the notation 384XY- 
B Z, where 384 refers to the number of grids in each spatial 
direction; X, Y = P, H3, H8 refers to physical or hyper- 
diffusion (and its power); Z=0, 1 refers to the strength of 
the external magnetic field. 

3.2. Parameter Space 

We require that, at the energy injection scale, B k^ L w 
Vtfc-UZ,. Therefore, we take B ~ 1 for most of the runs 
(see Table 1). 

For ordinary MHD turbulence, the dissipation cutoff 
(k c ~ I /Id) occurs at k ~ A/3. For viscosity-damped 
MHD turbulence, we require that k c <~ k^ (<~ 1/L) to 
maximize the dynamical range. For 384PH3-Bol, we use 
v <~ 0.015 (or, Reynolds number R <~ 100) to guarantee 
turbulence. In fact, the energy spectrum shows that the 
viscous cutoff occurs at k <~ 7 when we take v ~ 0.015. 
When we take v ~ 0.06 (e.g. run 256PH8-B 0.5), the vis- 
cous cutoff occurs right at the energy injection scale k^. 
In this case, we do not expect a turbulent velocity field. 

4. SPECTRA, ANISOTROPY, AND SCALE-DEPENDENT 
INTERMITTENCY 

4.1. Spectra and k\\ 

Figure 1 shows the energy spectra of ordinary and 
viscosity-damped turbulence. In both cases, energy is in- 
jected at k ~ 2.5. In ordinary turbulence, the injected en- 
ergy cascades down to the single scale at k <~ 100. In the 
viscosity-damped case, the kinetic energy and magnetic 
energy cascade together to the viscous damping scale at 
k ~ 7. Beyond k ~ 7, kinetic spectrum drops sharply due 
to viscous damping. In contrast the magnetic spectrum 
flattens out to a fc" 1 spectrum for 10 < k < 80. 

Figure 2 shows the parallel wave number k\\ as a func- 
tion of the total wavenumber. A method for calculating k\\ 
is described in CLV02b. The term B • VB describes mag- 
netic tension and is approximately equal to B^kubi for an 
isolated eddy in a uniform mean field Bo- Here fc|| oc 1/Zy 
and Z|| is the characteristic length scale parallel to B , 
which is known to be larger than the perpendicular length 
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1 10 100 1 10 100 

k k 

Fig. 1. — Spectra, (a: left) Ordinary MHD turbulence (384H8H8-Bol). The kinetic spectrum (E v (k)) is compatible with the Kolmogorov 
spectrum. However, the magnetic spectrum (Ef,(k)) is slightly shallower than Kolmogorov. (b: right) Viscosity-damped MHD turbulence 
(384PH3-Bol). The kinetic spectrum declines quickly after the viscous cutoff at k ~ 7. The magnetic spectrum follows a fc" 1 power-law after 
the cutoff, in agreement with theoretical expectations. The slope of E v (k) is somewhat steeper than the k~ 4 dependence that follows from 
theoretical considerations. Figure (b) is similar to that in CLV02b. 



scale. Cho & Vishniac (2000) and CLV02a argued that, in 
actual turbulence, eddies arc aligned with the local mean 
field B L . We can obtain the local frame representation of 
fcy , by considering an eddy lying in the local mean field Bl: 
Bl • Vb; w £?£fc|| &j. 3 The Fourier transform of this relation 

yields |Bl • Vb;| k ~ B^ku |b|k, where hatted variables are 
Fourier-transformed quantities. From this, we have 

1/2 



fell 



E 



fe<|k'|<fc+l 



|BL-Vb ; | k , 



Efc<|k'|<fc-| 



\*>\l 



(14) 



Figure 2(a) shows that, for ordinary turbulence, k\\ mea- 
sured by this method gives results consistent with the 

GS95 relations between k\\ and k\: h 



k 2 / 3 . This new 



method is complementary to the previous method utilizing 
structure functions (Cho & Vishniac 2000) for the study of 
anisotropy. Figure 2(b) shows that, for viscosity-damped 
turbulence, ku w constant. 

4.2. Scale- dependent intermittency 

The theoretical model in LVC03 predicts that intermit- 
tency is scale-dependent for viscosity-damped MHD tur- 
bulence. For ordinary MHD turbulence, no theory has 
addressed this issue. 

Figure 3 shows magnetic structures in a plane perpen- 
dicular to the mean field Bo- In Figures 3(a) and (b), we 
plot medium scale eddy shapes, using Fourier modes with 
10 < k < 20. (Remember that in the viscosity-damped 
turbulence viscous damping occurs around k <~ 7.) In Fig- 
ures 3(c) and (d), we plot small scale eddy shapes using 
Fourier modes with 40 < k < 80. In ordinary turbulence 
(Figures 3(a) and (c)), magnetic field structures are more 
or less smooth. In viscosity-damped turbulence (Figure 
3(b) and (d)), we can see that intermittency is more pro- 
nounced at small scales. It is obvious that intermittency 
is scale-dependent in viscosity-damped turbulence. 



Figure 4 shows scale-dependent intermittency more 
clearly. On the x-axis we plot the volume fraction and 
on the y-axis the fraction of the perturbed magnetic en- 
ergy contained in that volume. For example, in Figure 
4(b) we see that about 7% of the volume contains half the 
magnetic energy for modes with 40 < fc < 80. However, at 
larger scales, the same energy occupies larger volume (e.g. 
see the solid line). In ordinary strong MHD turbulence 
(Figure 4(a)), there are only slight differences between the 
curves corresponding to different scales. 

5. HIGH ORDER STATISTICS 

Spectra do not provide a full description of turbulence. 
Another useful source of information is the scaling of high 
order structure functions. The p-th order (longitudinal) 
velocity structure function SF P and scaling exponents ((p) 
arc defined as 

SF p (r) = (| [v(x + r) - v(x)] • r \ p ) oc r^ p \ (15) 

where the angle brackets denote averaging over x. In this 
section, we discuss the scaling of ((p) in ordinary and 
viscosity-damped turbulence. 

5.1. Ordinary turbulence 

The scaling relation suggested by She & Leveque (1994) 
contains three parameters (see Politano & Pouquet 1995; 
Muller & Biskamp 2000): g is related to the scaling 
vi <~ I 1 ! 9 , x related to the energy cascade rate tj 1 <~ l~ x , 
and C, the co-dimension of the dissipative structures: 

P, 



C(p) = J(i 



x) + C (l - (1 - x/C) p/9 ^j 



(16) 



For incompressible hydrodynamic turbulence, She & Lev- 
eque (1994) obtained 



e L (p)=p/9 + 2 l-(2/3) 



\P/S 



(17) 



3 There can be many ways to define the local mean field B^. In CLV02b and this paper, we obtain B^ for an eddy of (perpendicular) 
size I oc 1/k by eliminating modes whose perpendicular wavenumber is greater than k/2 and b; by eliminating modes whose perpendicular 
wavenumber is less than k/2. 
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Fig. 2. — Scaling of parallel wavenumber which is a measure of anisotropy. (a: left) Ordinary MHD. This result confirms the scale- 
dependent anisotropy prediction of Goldreich & Sridhar (GS95): fey ~ k 2 ^ 3 . (b: right) Viscosity-damped MHD. The parallel wavenumber fcy 
remains almost constant in the viscosity-damped regime (k > 10). Figure (b) is similar to one in CLV02b. 



using g = 3, x = 2/3, and C = 2, implying that dissipation 
happens over ID structures (e.g. vortices). 

In 3-dimensional MHD turbulence, there have been 
three recent developments for the scaling of the structure 
function exponents ((p). First, Miiller & Biskamp (2000) 
extended the She-Leveque model to incompressible MHD 
turbulence and obtained 

C MB (p)=p/9 + l-(l/3F/ 3 , (18) 
where they used g = 3, x = 2/3, and assumed that the dis- 
sipating structures are 2-dimensional (i.e C = 3 — 2 = 1). 
Their numerical calculations for decaying MHD turbu- 
lence without mean field show that scaling exponents for 
z ± = v ± b follow this relation very closely. Second, 
CLV02a studied MHD turbulence in a strongly magne- 
tized medium and showed that turbulent motions (i.e. 
velocity) in planes perpendicular to the local mean field 
directions follow the original She-Leveque scaling (equa- 
tion (17)). Physically this might mean that the dissipa- 
tive structures are again ID hydrodynamic-type vortices. 
However, more probably, this may just mean that the dis- 
sipation structures look like one-dimensional in the (per- 
pendicular) slices of 3D MHD turbulence. Third, Boldyrev 

(2002) assumed that in compressible turbulence dissipa- 
tion happens over 2D structures, i.e. C — 1 and eq. (18) 
for £. This scaling was supported by compressible MHD 
simulations (Padoan et al. 2003a). Padoan et al. (2003a) 
showed that the velocity field in highly super-sonic MHD 
turbulence follows the scaling given by eq. (18) scaling, 
while in subsonic MHD turbulence it follows the scaling in 
equation (17). 

Combined together these results look puzzling. One can 
see easily that there are at least two apparent contradic- 
tions. First, Cho et al.'s (CLV02a) result conflicts with 
Miiller & Biskamp's result. This discrepancy could be at- 
tributed to the different simulation parameters. For ex- 
ample, CLV02a used a strong mean field while MB00 used 
no mean field. However, according to Miiller & Biskamp 

(2003) , the difference in Bq does not resolve this problem. 
Second, Padoan et al. (2003a) is at odds with Miiller & 
Biskamp (2000): Padoan et al.'s results for low Mach num- 
ber do not converge to the incompressible case. The differ- 



ence in Bq looks marginal: Miiller & Biskamp (2000) used 
a zero mean magnetic field while Padoan et al. (2003a) 
used a weak mean magnetic field. 

These contradictions have led us to revisit the issue of 
intermittency in strong MHD turbulence. In Figure 5(a), 
we plot longitudinal velocity structure functions SF p . The 
calculations are done in the perpendicular planes in the lo- 
cal frame of reference, where the parallel axis is aligned 
with the local mean magnetic field (see CLV02a): 

SF local {r) = (| [v(x + r±) _ y(x)] . f± | P)) (1Q) 

where rj_ is perpendicular to local mean field and r± is 
the unit vector parallel to rj_. We show the second, third, 
fifth, and 10th order structure functions. We observe 
that the slopes are slightly shallower than those of She- 
Leveque model. Here the results are for longitudinal struc- 
ture functions, which mostly reflect Alfven mode statistics 
in the perpendicular planes. However, when we use nor- 
mal structure functions ((| [v(x + rj_) — v(x)] \ p }), which 
reflect both Alfven and pseudo- Alfven statistics, the C(p)' s 
show a larger deviation from the She-Leveque model. This 
may be understood if the pseudo- Alfven modes have dif- 
ferent scaling exponents. We can separate out the scaling 
of pseudo- Alfven modes when we use transverse structure 
functions in the perpendicular planes in local frame. We 
can also capture them using longitudinal structure func- 
tions in the parallel directions in local frame. However, we 
will not pursue this issue further here. 

In Figure 5(b), we plot the normalized differential slope, 
[d\n SFp/dlnr]/[d\n SFs/dlnr], for local frame velocity 
structure functions. In the actual calculations, we use 
[d In SF p /d In r](r) w \n[SF p (r + 4)/SF p (r - 4)]/ln[(r + 
4)/(r — 4)]. From the Figure, we obtain the ((p)'s by av- 
eraging the differential slope over r € [20,45]. In general, 
normalizing the differential slope using SF 3 gives better- 
defined scaling exponents, which was noted in Miiller & 
Biskamp (2003) and by Boldyrev (private communication). 

In Figure 5(c), we plot the ((p)'s obtained this way for 
velocity, magnetic field, and z ± . The scaling exponents for 
velocity show reasonable agreement with the She-Leveque 
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Fig. 3. — Magnetic structures in planes perpendicular to the mean magnetic field, (a) and (c): Ordinary MHD turbulence. The distribution 
of magnetic energy is more or less smooth. We do not observe scale-dependency, (b) and (d): Viscosity-damped MHD turbulence. In the Run 
384PH3-Bol, the viscosity-damped regime emerges below the dissipation scale at k c ~ 10 (see Figure 1(b)). Panel (b) is not much different 
from panel (a) because the corresponding scale is not far from the dissipation scale. In panel (d), the distribution of magnetic energy is very 
intermittent. Panels (b) and (d) illustrate that intermittency is more pronounced at small scales. 
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Fig. 4. — Fractional volume (X-axis) vs. fractional magnetic energy in the volume (Y-axis) (a: left) Ordinary MHD strong turbulence: 
intermittency is independent of scale, (b: right) Viscosity-damped MHD turbulence: smaller scales show a higher concentration of magnetic 
energy. 



model, in agreement with CLV02a. Those for the mag- 
netic field follow the Muller-Biskamp model for small p 
but level off and do not change much when p > 7. Those 
for Elsasser variables z ± lie between those for the velocity 
and magnetic fields. 

In Figure 5(d), we plot C(p)' s calculated in the global 
frame, in which coordinate axes are aligned with the usual 
Cartesian axes. The method we used for Figure 5(d) is 
similar to the methods used by Miiller & Biskamp (2000) 
and Padoan et al. (2003a), apart from the fact that our 
turbulence is strongly magnetized. It is worth noting that 
we get a good agreement with the Muller-Biskamp (2000) 
model for z ± variables. 

Our first result is that the scaling exponents obey the 
rule C maanetic (p) < ( z {p) < C veloclty {p). It matters 
whether one uses Elsasser variables, velocity or magnetic 
field to determine the dimension of the dissipation struc- 
tures. Our second result is that the dimension of the dissi- 
pation structures looks different when viewed in local and 
global frames of reference 4 . It looks as if one dimensional 
vortices merge into two-dimensional sheets when viewed 
from the global system of reference. These differences 
should make us wary of a naive association of the cor- 
responding parameters in cq.(16) with the dimensions of 
the dissipation structures in MHD turbulence. 

This result can eliminate the contradiction between 
Miiller & Biskamp (2000) and CLV02a: the former used 
the Elsasser variables in the global frame and the latter 
used velocity in the local frame for studying scaling expo- 
nents (see Figure 5(c) and (d)). Similarly, the conflict be- 
tween Miiller & Biskamp (2000) and Padoan et al. (2003a) 
can be partially relieved: the latter used velocity and, 
therefore, we expect that the scaling exponents in the lat- 
ter are larger than those in the former. However, it is not 
certain from this that the discrepancy between the Padoan 
et al. (2003a) scalings at low Mach number and those in 
Miiller & Biskamp (2000) in the incompressible case has 
been completely resolved. It is also unclear why the mag- 
netic and velocity fields have different scalings. The sense 
of the difference suggests that the magnetic field is signif- 



icantly more intermittent than the velocity field. Appar- 
ently this topic deserves further theoretical research. 

5.2. High order structure functions in decaying 
compressible MHD turbulence 

Above we speculated that there may not be a contradic- 
tion between the results by Miiller & Biskamp (2000) and 
Padoan et al. (2003a). To test this we performed a numer- 
ical simulation for decaying compressible MHD turbulence 
without mean magnetic field. The compressible MHD code 
is described in detail in Cho & Lazarian (2002). We use 
the same numerical scheme described above. The initial 
kinetic and magnetic energy spectra are 

E v {k) ^E b (k) «fc 2 e- fc2 / fc o, (20) 

where we take ko — 2. The initial phases of the individual 
Fourier components are random. The initial sonic Mach 
number (the ratio of the rms velocity to the sound speed) 
is ~ 0.7 and the density is constant. 

Initially only small- fc Fourier components near k — k 
are excited, but as the energy cascade begins to operate the 
large-fc Fourier components are excited. We measure high- 
order structure function statistics after the energy spectra 
develop inertial range (Figure 6(a)). Initially the velocity 
and magnetic fields have almost identical spectra (equa- 
tion (20)). It is interesting that the magnetic energy spec- 
trum is larger than its kinetic counterpart at later times 
(see Biskamp & Miiller 2000 for a similar result in the 
incompressible case). Figure 6(b) shows the raw velocity 
structure functions. Figure 6(c) shows the normalized £(p) 
(i.e.£(p)/£(3)), which we obtain by averaging the slopes 
between r = 17 and r = 34. Interestingly, it appears that 
the velocity closely follows She-Leveque scaling while the 
magnetic field follows Muller-Biskamp scaling. The veloc- 
ity scaling we find is indeed consistent with the result in 
Padoan et al. (2003a). However, since we do not resolve a 
sufficiently long inertial range (see the spectra), the results 
in Figure 6(c) should be regarded as very tentative. We 
will present a comprehensive study of higher order statis- 
tics elsewhere. 



4 Note that the calculation is done for the perpendicular planes in the local frame, while it is done for all directions in the global frame. 
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Fig. 5. — The intermittency of ordinary turbulence (384H8H8-Bq1). (a: upper-left) Velocity structure functions (multiplied by arbitrary 
constants) in planes perpendicular to the local mean magnetic fields, (b: upper-right) Differential slopes normalized by that of the third-order 
structure function: [d In SF P /d In r]/[d In SFg/d In r], where SF P is the p-th order longitudinal velocity structure function calculated in planes 
perpendicular to the local mean magnetic fields. The slope for SFi is the bottom curve and that for SF12 is the top curve, (c: lower-left) 
Normalized structure function exponents in perpendicular directions in the local frame. The velocity exponents show a scaling similar to the 
Shc-Leveque model. The magnetic field shows a different scaling, (d: lower-right) Normalized structure function exponents in the global 
frame. Note that the result for is very similar to the Miiller-Biskamp model. 
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Fig. 6. — Compressible decaying MHD turbulence (COM2). The mean magnetic field is zero: Bo=0. (a: left) Spectra, (b: middle) 
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5.3. Viscosity-damped turbulence 

The run 384PH3-Bol has a very limited inertial range 
for the magnetic field, and is not suitable for the study 
of scaling exponents. We use the run 256PH8-Bo0.5 
{y = 0.06 and r\ = hyper — diffusion) instead. In this 
run, the mean field strength is reduced to 0.5 because 
the rms velocity is ~0.4 due to strong viscous damping. 
The energy spectra is plotted in Figure 7. This run also 
clearly shows scale dependent intermittency. The paral- 
lel wavenumber k\\ shows a very slow increase, a factor of 
~ 2 from k = 5 to k = 70. Viscous damping occurs right 
at the energy injection scale (Figure 7(a)). The magnetic 
spectrum shows a slope roughly consistent with —1. Here 
we consider only the magnetic scaling exponents. 

Figure 7(b) shows the structure functions for magnetic 
field. The calculation is done in perpendicular planes in 
local frame. Figure 7(c) shows that the scaling exponents 
becomes negative when p is larger than 8. 

In the original She-Leveque model, they first obtained 
the scaling relation for high order statistics of energy dis- 
sipation averaged over a ball of size I, ef. < ef >^ l Tp with 
t sl = _2p/3 + 2[l- (2/3) p ]. Then, using the Kolmogorov 

refined similarity hypothesis, e/ vf/l or vi w i 1 / 3 ^ 3 , 
they obtained 

( SL (p)=p/3 + r p/3 . (21) 

Note that t p < and, hence, ( SL (p) < p/3. In viscosity- 
damped MHD turbulence, we cannot use the Kolmogorov 
refined similarity hypothesis directly. Instead, from ei w 
bf/l° or bi w e 1 / 2 , we expect 

C damped (p) = r p/2 . (22) 

As in She-Leveque (1994; see also Politano & Pouquet 
1995), we may write 

r p = -xp + C[l-{l-x/C) p ], (23) 
where x and C are defined as in equation (16). When we 
substitute x = 0, we get t p = and Ramped ^ = Q How _ 
ever, the 2-parameter She-Leveque model in equation (16) 
may not be applicable to the case of x = 0. 

Despite the evident problems with extending the She- 
Leveque model into the viscosity-damped regime, a crude 
estimate for the £(p) can be obtained from the theoretical 
model discussed in §2. Given the extreme intermittency of 
the magnetic field we expect that equation (15) will give 

SF p (l) ~ bffr ~ l 1 -^ 2 , (24) 

or 

C(P) = 1-|- (25) 

However, the rise in magnetic field strength at small length 
scales introduces a competing contribution from the resis- 
tive scale, since for p > 2 the rise in local magnetic field 
strength overwhelms the decrease in the filling factor as 
a function of scale. This small scale contribution is not a 
function of I and when it dominates we expect ((p) <~ 0. 
Putting all this together, we see that our simple model 
for the viscosity-damped regime predicts that ((1) will lie 
between and 0.5, ((2) <~ 0, and when p > 2 we expect 
that ((p) will become negative and then asymptote to 
at large p. 

Comparing this prediction to Figure (7) we see that this 
naive prediction enjoys some qualitative success, but is cer- 
tainly not exactly correct. Part of the problem may lie in 



the limited resolution of our numerical simulation. For ex- 
ample, the expectation that the magnetic energy spectrum 
follows a fc _1 scaling is fairly compelling, but in that case 
we would expect £(2) = 0. The scaling in the simulation 
is slightly steeper (Figure 7(a)). The failure of the higher 
order exponents to match the theoretical model may be 
part of the same problem, but it is also likely to reflect 
the crude nature of the model. Kolmogorov turbulence 
theory gives a passable match to C(2), but fails noticeably 
at higher p, and our model for the viscosity-damped MHD 
regime is not particularly more sophisticated. 

6. TESTING AND DISCUSSION 

6.1. Viscosity- damped turbulence with physical diffusion 

How real are the structures that we observe in the new 
regime of MHD turbulence? Could they be a numerical 
artifact? 

We have used a third or higher order hyper-diffusion 
term to minimize the effects of magnetic diffusion. How- 
ever, in general, high order hyper-diffusion suffers from a 
bottle-neck effect, which is characterized by a flattening of 
the energy spectrum near the dissipation scale. Therefore 
it is necessary to check if the tail of the magnetic fluctua- 
tions is a real physical effect and not due to a bottle-neck. 
Here, we present a run with a physical magnetic diffusion 
term (Run 256PP-B 1; v = 0.015, r, = 0.001). 

In Figure 8(a), we see that, compared with its kinetic 
counterpart, the magnetic spectrum has more power at 
k > 10. Although we cannot say much about the spec- 
trum of the magnetic fluctuations at scales smaller than 
the viscous cutoff (k c ~ 7), the existence of the new regime 
is evident in Figure 8(a). A visual inspection of the mag- 
netic structures in Figure 9(a) exhibit intermittent struc- 
tures similar to those seen in the case of hyper-diffusion. 
We see that the bottle-neck effect is not responsible for 
intermittent magnetic structures at small scales. 

Further testing of the viscous damped regime was per- 
formed in Cho & Lazarian (2003a, b). There, using the 
compressible code mentioned earlier, we calculated not 
only velocity and magnetic field strength, but also den- 
sity fluctuations that arise from MHD turbulence in the 
viscosity-damped regime. We present our results in Figure 
8(b), side by side with a plot of the incompressible simu- 
lation results. The similarity between the two are vivid. 
It is also evident that the shallow spectrum of magnetic 
fluctuations results in a shallow spectrum of density fluc- 
tuations. Theoretically we expect the density spectrum to 
follow a fc _1 law. Unfortunately, our low resolution does 
not allow us to test this prediction at the moment. 

Figure 9(c) shows that the density also has intermittent 
structures below the viscous cutoff. The density fluctu- 
ations may be (anti-)correlated with magnetic structures 
(see Figure 9(b)). However, we will postpone further dis- 
cussion of this point until higher resolution runs become 
available. 

6.2. Comparison with Observations 

In this paper, we have considered properties of ordinary 
and viscosity-damped MHD turbulence. It is important to 
test our predictions with observations. 
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Fig. 8. — Spectra of runs without hyper-diffusion. We use the same physical viscosity (y = 0.015) for velocity, (a: left) From an incom- 
pressible run with a physical magnetic diffusion (256PP-Bol). We use a small physical magnetic diffusion coefficient (rj = 0.001). The kinetic 
spectrum declines quickly due to the large viscosity. The magnetic spectrum shows structures at small scales (i.e. k > 10). (b: right) From 
a compressible MHD run (COM1). We use numerical diffusion. The rms Mach number is ~0.5. Due to the large viscosity (v = 0.015), the 
kinetic spectrum drops quickly for k > 7. The magnetic spectrum shows structures at small scales (i.e. k > 10). Note that the shape of the 
density spectrum (E p ) is similar to the magnetic spectrum. Figure (b) is from Cho & Lazarian (2003a,b). 




Fig. 9. — Intermittent small scale structures in a plane perpendicular to the mean magnetic field, (a: left) From an incompressible run 
with a physical magnetic diffusion (256PP-B()1). Small scale (k > 20) magnetic field structures show intermittency. (b: middle) From a 
compressible run (COM1). Small scale (k > 20) magnetic field structures, (c: right) From COM1. Small scale [k > 20) density structures. 
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Indirect Evidence — Existence of density structures at 
scales smaller than a fraction of a parsec is particularly im- 
portant for the small scale density structures in the ISM. 
We speculate that the power-law density tail might have 
some relation to the tiny-scale atomic structures (TSAS). 
Heilcs (1997) introduced the term TSAS for the mysterious 
H I absorbing structures on the scale from thousands to 
tens of AU, discovered by Deiter, Welch & Romncy (1976). 
Analogs are observed in Nal and Call (Meyer & Blades 
1996; Faison & Goss 2001; Andrews, Meyer & Lauroesch 
2001) and in molecular gas (Marscher, Moore & Bania 
1993). Recently Deshpande, Dwarakanath & Goss (2000) 
analyzed channel maps of opacity fluctuations toward Cas 
A and Cygnus A. They found that the amplitudes of den- 
sity fluctuations at scales less than 0.1 pc are far larger 
than expected from extrapolation from larger scales, con- 
sistent with the existence of TSAS. 

Measuring velocity and density spectra — Velocity chan- 
nel analysis (VCA; Lazarian & Pogosyan 2000) is a recent 
technique that can extract velocity and density spectra. 
The technique predicts that the power spectra of spectral 
line intensity vary when the thickness of velocity chan- 
nels changes 5 . The compressible simulation in the previous 
subsection suggests that density may have a k~ x spectrum 
below the viscous cutoff. At these scales, velocity effects 
should be suppressed. Deshpande et al. (2000) used ab- 
sorption measurements to test the HI structure on subpar- 
sec scales and did not see the variations of the channel map 
spectra as they added their channels together. This meant 
that the velocity effects were negligible at the scales they 
studied. According to the VCA this meant that the shal- 
low spectrum of density k~ a (a < 1) that they measured 
corresponds to real density fluctuations. This spectrum 
may arise from the new regime of MHD turbulence. Fur- 
ther tests, including those which use different techniques 
are necessary. 

Information on the velocity statistics can also be ob- 
tained using centroids of velocity (see Munch 1958). Mod- 
ified velocity centroids (MVCs) that are not sensitive to 
density fluctuations and can therefore be used for studies 
of velocity statistics were proposed in Lazarian & Esquivcl 
(2003). They argued that combining the VCA and MVCs 
it is possible to get reliable measures of the underlying 
velocity spectra. It will be important to analyze the Desh- 
pande et al. data using MVCs. 

Studies of anisotropy Testing scale-dependent 
anisotropy is a challenging problem, because scale- 
dependent anisotropy is averaged away when we observe 
turbulence from outside. The scale-dependent anisotropy 
is revealed only in a local frame of reference whose paral- 
lel axis is parallel to the local mean magnetic field (Cho & 
Vishniac 2000; CLV02a). Nevertheless, we can still study 



global anisotropy with 3D maps using velocity centroids 
(Lazarian, Pogosyan & Esquivel 2002) or 2D maps ob- 
tained from the Position-Position- Velocity data (Esquivel 
et al. 2003). 

Higher order statistics — Scale-dependent intermittency 
is hard to study using observational data. Usually the 
noise increases rapidly as we go to higher order statistics. 
However, high order density structure function scalings 
have been used to study interstellar turbulence by Padoan 
et al. (2003b) and Padoan, Cambresy, & Langer (2003). 
Further studies of the problem are obviously required. 

7. SUMMARY 

We have compared ordinary MHD turbulence and 
viscosity-damped MHD turbulence. They are different in 
many ways - in their spectra, in the degree of anisotropy 
and its scaling with length, and in the nature of intermit- 
tency in the turbulent cascade. 

1. The spectra of ordinary strong MHD turbulence 
are compatible with Kolmogorov's /c~ 5 / 3 spectrum. 
The magnetic spectra in viscosity-damped turbu- 
lence show a flatter spectrum: E(k) ~ fc _1 . 

2. Eddies in ordinary MHD turbulence show scale- 
dependent anisotropy consistent with the Goldre- 
ich & Sridhar model. Eddies in viscosity-damped 
MHD turbulence show extremely anisotropic struc- 
tures - their parallel size does not decrease signifi- 
cantly with the perpendicular eddy size. 

3. Intermittency of magnetic field strength in ordinary 
turbulence is not scale-dependent. Intermittency 
in viscosity-damped turbulence is scale-dependent: 
smaller scales are more intermittent. 

4. The velocity and magnetic fields show different scal- 
ing exponents for high order structure functions. 

5. Scaling exponents for the magnetic field in 
viscosity-damped MHD turbulence show very little 
change as order of the structure changes. They are 
generally small and change from positive at p = 1 
to negative values at moderate p. 
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Table 1 
Simulations 



Run a 


N 3 


V 


V 


Bo 


Note 


384H8H8-B 1 


384 3 


hyper 


hyper (V 16 ) 


1 


incompressible 


384PH3-B 1 


384 3 


.015 


hyper (V 6 ) 


1 


incompressible 


256PH8-Bo0.5 


256 3 


.06 


hyper (V 16 ) 


.5 


incompressible 


256PP-B 1 


256 3 


.015 


.001 


1 


incompressible 


COM1 


216 3 


.015 


numerical 


1 


compressible 


COM2 


256 3 


numerical 


numerical 





compressible 



a We use the notation 384XY-B Z (or 256XY-S0Z), where 384 (or 256) 
refers to the number of grid points in each spatial direction; X, Y = P, H3, 
H8 refers to physical or hyper-diffusion (and its power); Z=0.5, 1 refers to 
the strength of the external magnetic field. 



